Simple and faithful nonlinear field equations for aligning self-propelled rods 
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We derive a set of minimal yet complete nonlinear field equations describing the collective proper- 
ties of self-propelled rods from a simple microscopic starting point, the Vicsek model with nematic 
alignment. Analysis of their linear and nonlinear dynamics shows good agreement with the original 
microscopic model. In particular, we derive an explicit expression for the fronts forming density- 
segregated, banded solutions, allowing us to develop a more complete analytic picture of the problem 
at the nonlinear level. 
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Collective motion is a central theme in the rapidly 
growing field of active matter studies which loosely 
groups together all situations where energy is spent lo- 
cally to produce coherent motion [l| . In spite of the emer- 
gence of better-controlled, larger-scale experiments [Il- 
ly] , our understanding of collective motion mostly comes 
from the study of mathematical models, and particularly 
of models of "dry" active matter systems, where the fluid 
which surrounds the moving objects can be neglected. 

Microscopic models then usually consist of interacting 
self-propelled particles, as in the Vicsek model Q, where 
constant-speed point particles "ferromagnetically" align 
their velocities with that of local neighbors. The study 
of these models, together with some more theoretical ap- 
proaches, revealed a wealth of phenomena such as true 
long-range orientational order in two dimensions, sponta- 
neous segregation of dense/ordered regions, anomalously 
strong number fluctuations, etc. [8l-ll0l|. 

These results have given rise to an emerging picture of 
universality classes, typically depending on the symme- 
tries involved, which one would ideally characterize by 
some coarse-grained field equations. Different routes can 
be followed to obtain such equations: they can be written 
a priori, putting down all terms allowed by symmetries 
up to some arbitrary order in gradients once the hydro- 
dynamic fields have been identified. One can also de- 
rive continuous descriptions from some microscopic start- 
ing point under more or less controlled and constrain- 
ing assumptions, yielding more or less complete, well- 
behaved equations. There is nevertheless shared belief, 
based mostly on renormalization-group approaches, that 
in each case there exists a set of minimal equations ac- 
counting for all large-scale physics. 

For polar particles aligning ferromagnetically (as in the 
Vicsek model), there is now near-consensus about this 
minimal set of equations: the phenomenological theory, 
initially proposed by Toner and Tu fill ] is essentially cor- 



rect if one takes into account the dependences of its coeffi- 
cients on density and parameters initially overlooked but 
later derived from microscopies in [ijl, EH . It has been 
shown to reproduce many of the phenomena observed in 
microscopic models [IH , although a complete study of its 
nonlinear solutions and dynamics is still lacking. 

For the other important universality class of polar par- 
ticles aligning nematically — e.g. self-propelled rods — , 
the situation is less satisfactory: Baskaran and Marchetti 
[IH first derived rather lengthy yet mostly linear hydro- 
dynamic equations for hard rods interacting via excluded 
volume, showing in particular the presence of global ne- 
matic, not polar, order, in agreement with microscopic 
observations [lj|. Very recently [It}, they added some 
nonlinear terms and performed a linear stability analy- 
sis of the homogeneous ordered state within an arbitrary 
choice of parameters. These even longer equations do not 
benefit from the same consensus as the Toner- Tu theory. 

In this Letter, we derive a set of minimal yet complete 
nonlinear field equations describing the collective prop- 
erties of self-propelled rods from a simple microscopic 
starting point, the Vicsek model with nematic alignment 
studied in [l6(. We use a "Boltzmann-Ginzburg-Landau" 
approach, a controlled expansion scheme 18| which is a 
refined version of that used in [l2| ■ Analysis of the linear 
and nonlinear dynamics of these equations shows good 
agreement with the original microscopic model. In par- 
ticular, we derive explicit expressions for the fronts form- 
ing density-segregated, banded solutions, allowing us to 
develop a more complete analytic picture of the problem 
at the nonlinear level. 

The Vicsek model with nematic alignment [3] consists 
of point particles moving off-lattice at constant speed Un- 
Orientations and positions are updated following (here in 
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where e(0) is the unit vector along 6, the sum is taken 
over particles k within distance do of particle j (includ- 
ing j itself), and 77 is a white noise with zero average and 
variance a 2 . Like all Vicsek-style models, it shows orien- 
tational order at large-enough global density p g an d/or 
small-enough noise strength a. It was shown in [16[ that 
the order is nematic and that both the ordered and dis- 
ordered phases are subdivided in two: The homogeneous 
nematic phase observed at low noise is replaced at larger 
a values by a segregated phase where a dense, ordered 
band occupying a fraction of space coexists with a dis- 
ordered, dilute, gas. The transition to disorder is given 
by the onset of a long- wavelength instability of this band 
observed upon increasing a further. Finally, this chaotic 
regime where dense ordered bands constantly form, elon- 
gate, meander, and disappear over very long timcscales is 
finally replaced by a "microscopically-disordered" phase 
at large noise intensities. 

Following [l2j , we write, in a dilute limit where only bi- 
nary interactions are considered and assuming that orien- 
tations are decorrelated between them ( "molecular chaos 
hypothesis"), a Boltzmann equation governing the evo- 
lution of the one-particle distribution f(r,0,t): 

d t f(r,e,t) + v e{6) ■Wf(r,6,t) = 1^/] + I col [f] (3) 

with the angular diffusion and collision integrals 
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where Pa{rj) is the microscopic noise distribution, <5 27r is 
a generalized Dirac delta imposing that the argument is 
equal to zero modulo 2ir, K (61,62) = 2doVo\e(6i) — e(0 2 )| 
is the collision kernel for dilute gases [IH, and * (61,62) = 
5 (0i + 2 )+f(.ff[cos(0 1 -0 2 )]-l) for -f < 02-01 < % 
(with H(x) the Heaviside step function) codes for the 
nematic alignment. Rescaling of time, space and density 
allows us to set the "collision surface" S = 2dov /X = 1 
and vq = 1 below, without loss of generality. 

Next, the distribution function is expanded in Fourier 
series of the angle: f(r,6,t) = ± fk(r,t)e~ lke , 
with fk = f *_ k and |/fe| < fo- The zero mode is nothing 
but the local density, while fi and f 2 give access to the 
polar and nematic order parameter fields P and Q: 
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As a matter of fact, it is convenient to use fi and /2, 
together with the "complex" operators V = d x + id y and 
V* = d x — idy. The continuity equation governing p is 
given by integrating the Boltzmann equation over angles: 

d i/9 + Rc(V*/i) =0. (6) 

In Fourier space, the Boltzmann equation ([3]) yields an 
infinite hierarchy of equations: 
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To truncate and close this hierarchy, we adopt the fol- 
lowing scaling structure, valid near onset of nematic or- 
der, assuming, in a Ginzburg-Landau-likc approach 18J, 
small and slow variations of the density and of the polar 
and nematic fields: 



p-po^e, {fok-l, / 2 fc}fc>i ~ e A 
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Note that the scaling of space and time is in line with 
the propagative structure of our system, as seen in the 
continuity equation ^ , which contains no diffusion term. 

The first non-trivial order yielding well-behaved equa- 
tions is e 3 : keeping only terms up to this order, equations 
for /fe>4 identically vanish, while those for fa and fi pro- 
vide expressions of these quantities in terms of p, fi , and 
/2, which allows us to write the closed equations: 



dtfi = ~(Vp + vVa) + ^/ 2 *V/2 



>-/3|i2| 2 )/i + C/i7 2 



(10) 



dth 



-2 v/i + -rV v*/ 2 - 5 A*v/ a - I V(/i/ 2 ) 



+ (M-ei/ 2 | 2 )/ 2 +^/ 1 2 + r|/ 1 | 2 / 2 



(11) 



where all coefficients depend only on the noise strength 
a (via the P^ coefficients) and the local density p: 
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FIG. 1: (color online) (a) Linear stability of homogeneous so- 
lutions in the (po, o) plane (plotted as a function of a 2 to en- 
hance clarity). The line p = — }^(P2-i) — given by u = 0, 

is the basic instability line defining pt(cr) or a t (po) : above 
it, the disordered homogeneous solution is linearly stable; be- 
low, it becomes unstable and the ordered solution fa = ^/ p/£ 
exists. This solution is unstable between the ot and the a s 
lines. It is linearly stable between a s , and a u , which marks 
the border of a region where q — is the most unstable mode. 
The color scale codes for the angle between the most unstable 
wavevector and the direction of nematic order, (b) Largest 
eigenvalue s+ (when positive) as a function of a 2 for po = 1. 



Below, we specialize for convenience to the Gaussian 
noise distribution P a {rf) = a ^ 2v ex p[~ f° r which 

A; = exp[— ifc 2 cr 2 ] [22||. A few remarks are then in order. 
First, p, can change sign and become positive for large 
enough p, while a is always positive: The homogeneous 
disordered state (/i=/2 = 0) undergoes an instability to 
nematic order when p, = 0, defining the basic transition 
line cr t (p ) in the (po,c) plane (Fig. QJi). Next, £ being 
positive in the p > region where the disordered so- 
lution is unstable, Eqs. (|10llll|) possess a homogeneous 
nematically-ordered solution (/i,/2) = (0, y/^TO (as- 
suming the order along x, so that fa is real positive) 

& 

In the following, we restrict ourselves to slightly mod- 
ified coefficients by suppressing the local density depen- 
dence (replacing p by po) of v and £, and thus of all other 
coefficients but the crucial linear ones p and a (note that 
this corresponds to retaining only terms up to order e 3 
after expanding coefficients in Sp). This expansion does 
not change any of the main features of our equations, but 
allows us to find exact front solutions (see below). 

We have studied the linear stability of the homoge- 
neous nematic solution with respect to perturbations of 
arbitrary wavevector in the full (po,o~) parameter plane 
(Fig. [IJ. Similarly to the polar case with ferromag- 
netic alignment 12-14, [2^, this solution is unstable to 
long wavelengths in a region bordering the basic transi- 
tion line. The most unstable modes in this region are 
roughly — but not exactly — transversal to the order of 
the solution [HI . The homogeneous nematic solution be- 
comes linearly stable deeper in the ordered phase (line 
<7 S in Fig. OJ,), but its stability domain is limited by an- 
other instability region where q — is the most unsta- 




FIG. 2: (color online) Numerically obtained density- 
segregated solutions, (a) density and fa profiles of a station- 
ary banded solution (fa = throughout). The fronts linking 
the disordered and ordered domains can be perfectly fitted to 
hyperbolic tangents (not shown). (cr = 0.26, po = 1,L = 200) 
(b,d,e) chaotic band regime: snapshots of (respectively) p, 
|/i I, and |/2 1. (a = 0.2826, po = 1,L = 200) 



ble mode (line er u , which can be shown to be given by 
a + C\/ m/£ — /5/V£ = 0)- This strong instability, which 
occurs at large densities and/or weak noise, may be an 
artifact introduced by our truncation [26l |. 

To go beyond the linear stability analysis of spatially 
homogeneous solutions, we performed numerical integra- 
tions of Eqs. (|10|llj) in rectangular domains with periodic 
boundary conditions of typical linear sizes 50-200 pBT ]. 
For parameter values in the region of transversal insta- 
bility of the nematic homogeneous solution, we observe, 
after transients and from almost every initial condition, 
stationary asymptotic solutions in which nematic order 
is confined to and oriented along a dense band with lo- 
cal density pband amidst a homogeneous disordered "gas" 
with /9 gas such that p ban d > Ps > Pt > Pgas (Fig. [2ji), 
where p s (cr) is given by inverting a s (po)- Varying system 
size and using various domain aspect ratios, we find most 
often a single band oriented along the shortest dimension 
of the domain, which occupies a size-independent fraction 
fl of space. All these observations are in agreement with 
the behavior of the original microscopic model [l6| . 

Band solutions are also present beyond the parameter 
domain where the homogeneous ordered state is linearly 
unstable. Starting from such a solution found at other 
parameter values — or from sufficiently inhomogeneous 
initial conditions — , we find band solutions both for po 
values larger than p s , where they coexist with the homo- 
geneous ordered phase, and below p t , where the disor- 
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FIG. 3: (color online) Analytic band solutions for the slightly 
simplified system (see text), (a) (po,f) parameter plane with 
basic instability line cr t , stability limit of homogeneous or- 
dered phase a s , and limits of existence of band solutions cr m i n 
and Umax- The short-dashed blue and red lines show the /9 gas 
and pband density values of the band solutions for po — 1 as a 
function of a over their existence range [er m m, er max ], indicated 
by the thin horizontal dashed violet lines, (b) variation with 
a of $7, the fraction of space occupied by the ordered part of 
the band solution, for po = 1. 



dered homogeneous solution is linearly stable. Working 
at fixed po varying a for clarity, we thus find bands in 
a [cmin, "max] interval larger than the linear- instability 
interval [<r 8 , <7t] (Fig. • Along it, the fraction occupied 
by the ordered band decreases from SI < 1 near <7 m j n 
to SI > near <r max . Furthermore, within a small layer 
"max < cr < a c , the bands are unstable, giving rise to 
a chaotic regime where they twist, elongate, break, and 
form again, in a manner strikingly similar to observations 
made in the original microscopic model (Fig. [2b, 123). 

Thus the region of linear instability of the homoge- 
neous ordered solution does not correspond to the exis- 
tence (and stability) domain of band solutions, which is 
wider. In the original microscopic model, with its built- 
in fluctuations, coexistence of band solutions and homo- 
geneous order has not been reported, but the homoge- 
neous solution was found metastable near the threshold 
of emergence of bands where these appear "suddenly" 
fli| . At the other end of the band existence region, no 
coexistence was reported between band solutions (chaotic 
or not) and the homogeneous disordered state, suggesting 
that the latter is always driven to the former by intrinsic 
fluctuations. Although this calls for revisiting the mi- 
croscopic model, this suggests that its thresholds do not 
correspond to those given by the linear stability of the 
homogeneous solutions of our continuous equations. 

We now show that many of the features of band so- 
lutions can be found analytically. Suppose that, as ob- 
served, fi — for band solutions and that f% is real and 
positive (i.e., nematic order is along x), and depends only 
on y. For a stationary solution, Eq. (|TU|) then yields, after 
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where p is a constant. This allows us to write Eq. (jlll) . 
again looking for stationary solutions, in terms of f% only: 
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where we have rewritten p = p'(p — p t ), with p! inde- 
pendent of p. Eq. (|14[) can be interpreted as an ordinary 
differential equation describing the motion in fictitious 
time y of a particle of mass j in a potential with one 
maximum at f 2 = and another at some value fi (close 
enough to threshold). A band solution is found for a tra- 
jectory starting at and returning to the same fixed point 
fi = (homoclinic orbit). Direct integration of Eq. (fT4"]) 
yields, under the condition lim^-^ioo f2(y) = 0, the fol- 
lowing solution 
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where a = \Jl + 9&(/5 — p t )/2p' and b = £ — p'j/2. 
The value of p can be obtained from the condition 
J L p (y) dy = Lpo, where L is the size of the box. We can 
neglect the exponentially decreasing tails in the integral 
and solve the equation J_ (p (y) — p) dy = L (po — p). 
Under the assumption L — > oo we obtain 
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where K\ and Ki are positive quantities depending on 
<7 and po whose expression we omit for compactness. 
Substituting this value in the expression of a gives us 
a = K\e~ K ^ L l 2 , yielding a width of the band propor- 
tional to L, in agreement with observations on the mi- 
croscopic model. As L — > oo, the value of p converges to 
the asymptotic value p gas - 

To determine the surface fraction SI occupied by the 
ordered band, we use the relation S7 (pband ~ /°gas) + 
Pgas = Po- Substituting the value of pband obtained from 
Eqs. (JTB"]) and ([To]) at y = 0, we find for L — > oo 
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The condition < SI < 1 yields the lower limit a m - ln and 
the upper limit cr max of existence of the band solution. 
All these results are presented in Fig. [3] 

To summarize, using a "Boltzmann-Ginzburg-Landau" 
controlled expansion scheme, we derived a set of mini- 
mal yet complete nonlinear field equations from the Vic- 
sek model with nematic alignment studied in (l6j . This 
simple setting allowed for a comprehensive analysis of 
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the linear and nonlinear dynamics of the field equations 
obtained because our approach automatically yields a 
"meaningful manifold" parameterized by global density 
and noise strength in the high-dimensional space spanned 
by all coefficients of the continuous equations. Excellent 
semi-quantitative agreement was found with the simu- 
lations of the original microscopic model. An analytic 
picture of the nonlinear banded solutions present was 
obtained. Their existence domain was found different 
from the region of linear instability of the homogeneous 
ordered phase, stressing the importance of a nonlinear 
analysis. More work, beyond the scope of this Letter, is 
needed to obtain a better understanding of the chaotic 
regimes observed. To this aim, we plan to study the lin- 
ear stability of the band solutions in two dimensions. 

Our equations (jlOlllip are simpler than those written 
by Baskaran and Marchetti in 15|, [l7[ , not only because 
our microscopic starting point does not include positional 
diffusion of particles. The method used there seems in- 
trinsically different, yielding more terms, many with a 
different structure from ours, while some of our nonlin- 
ear ones do not appear. It is not known yet whether 
the equations of [ljj can also account for the nonlincar 
phenomena described above. Future work should explore 
this point in some detail. 

Finally, recent experiments have revealed that micro- 
tubules displaced by a "carpet" of dynein motors grafted 
to a substrate self-organize their collective motion into 
large vortical structures (fjj. It was shown that a Vicsek 
model with nematic alignment, but one in which the mi- 
croscopic noise is colored, accounts quantitatively for the 
observed phenomena. Extending the approach followed 
here to this case is the subject of ongoing work. 
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